{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "IN_COLAB = 'google.colab' in sys.modules\n",
    "if IN_COLAB:\n",
    "    !git clone https://github.com/cs357/demos-cs357.git\n",
    "    !mv demos-cs357/figures/ .\n",
    "    !mv demos-cs357/additional_files/ ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hilbert Matrix - Condition Number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy.linalg as la\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ndim = np.array([2,3,8,11,14])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's perform linear solves for matrices with increasing size \"n\", for a problem in which we know what the solution would be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for nd in ndim:\n",
    "    ## This is the vector 'x' that we want to obtain (the exact one)\n",
    "    x = np.ones(nd)\n",
    "    ## Create a matrix with random values between 0 and 1\n",
    "    A = np.random.rand(nd,nd)\n",
    "    ## We compute the matrix-vector multiplication \n",
    "    ## to find the right-hand side b\n",
    "    b = A @ x\n",
    "    ## We now use the linear algebra pack to compute Ax = b and solve for x\n",
    "    x_solve = la.solve(A,b)\n",
    "    ## What do we expect? \n",
    "    print(\"------ N =\", nd, \"----------\")\n",
    "    error = x_solve-x\n",
    "    print(\"Norm of error = \", la.norm(error,2)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(x_solve)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will perform the same computation, but for a special matrix, known as the Hilbert matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Hilbert(n):\n",
    "    \n",
    "    H = np.zeros((n, n))    \n",
    "    for i in range(n):        \n",
    "        for j in range(n):        \n",
    "            H[i,j] = 1.0/(j+i+1)    \n",
    "    return H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for nd in ndim:\n",
    "    ## This is the vector 'x' that we want to obtain (the exact one)\n",
    "    x = np.ones(nd)\n",
    "    ## Create the Hilbert matrix\n",
    "    A = Hilbert(nd)\n",
    "    ## We compute the matrix-vector multiplication \n",
    "    ## to find the right-hand side b\n",
    "    b = A @ x\n",
    "    \n",
    "    ## We now use the linear algebra pack to compute Ax = b and solve for x\n",
    "    x_solve = la.solve(A,b)\n",
    "    ## What do we expect? \n",
    "    print(\"------ N =\", nd, \"----------\")\n",
    "    error = x_solve-x\n",
    "    print(\"Norm of error = \", la.norm(error,2)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(x_solve)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What went wrong?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Condition number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution to this linear system is extremely sensitive to small changes in the matrix entries and the right-hand side entries. What is the condition number of the Hilbert matrix?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for nd in ndim:\n",
    "    ## This is the vector 'x' that we want to obtain (the exact one)\n",
    "    x = np.ones(nd)\n",
    "    ## Create the Hilbert matrix\n",
    "    A = Hilbert(nd)\n",
    "    ## We compute the matrix-vector multiplication \n",
    "    ## to find the right-hand side b\n",
    "    b = A @ x\n",
    "    ## We now use the linear algebra pack to compute Ax = b and solve for x\n",
    "    x_solve = la.solve(A,b)\n",
    "    ## What do we expect? \n",
    "    print(\"------ N =\", nd, \"----------\")\n",
    "    error = x_solve-x\n",
    "    print(\"Norm of error = \", la.norm(error,2)) \n",
    "    print(\"Condition number = \", la.cond(A,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Residual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for nd in ndim:\n",
    "    ## This is the vector 'x' that we want to obtain (the exact one)\n",
    "    x = np.ones(nd)\n",
    "    ## Create the Hilbert matrix\n",
    "    A = Hilbert(nd)\n",
    "    ## We compute the matrix-vector multiplication \n",
    "    ## to find the right-hand side b\n",
    "    b = A @ x\n",
    "    ## We now use the linear algebra pack to compute Ax = b and solve for x\n",
    "    x_solve = la.solve(A,b)\n",
    "    ## What do we expect? \n",
    "    print(\"------ N =\", nd, \"----------\")\n",
    "    error = x_solve-x\n",
    "    residual = A@x_solve - b\n",
    "    print(\"Error norm = \", la.norm(error,2)) \n",
    "    print(\"Residual norm = \", la.norm(residual,2)) \n",
    "    print(\"Condition number = \", la.cond(A,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rule of thumb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for nd in ndim:\n",
    "    ## This is the vector 'x' that we want to obtain (the exact one)\n",
    "    x = np.ones(nd)\n",
    "    ## Create the Hilbert matrix\n",
    "    A = Hilbert(nd)\n",
    "    ## We compute the matrix-vector multiplication \n",
    "    ## to find the right-hand side b\n",
    "    b = A @ x\n",
    "    ## We now use the linear algebra pack to compute Ax = b and solve for x\n",
    "    x_solve = la.solve(A,b)\n",
    "    ## What do we expect? \n",
    "    print(\"------ N =\", nd, \"----------\")\n",
    "    error = x_solve-x\n",
    "    residual = A@x_solve - b\n",
    "    print(\"Error norm = \", la.norm(error,2)) \n",
    "    print(\"Residual norm = \", la.norm(residual,2)) \n",
    "    print(\"|dx| < \", la.norm(x)*la.cond(A,2)*10**(-16))\n",
    "    print(\"Condition number = \", la.cond(A,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rule of Thumb on Conditioning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# print(plt.style.available) # uncomment to print all styles\n",
    "import seaborn as sns\n",
    "sns.set(font_scale=2)\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "plt.rcParams['figure.figsize'] = (8,6.0)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's make a matrix\n",
    "\n",
    "Make the second column nearly linearly indepent to the first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = 10\n",
    "A = np.random.rand(n,n)\n",
    "\n",
    "delta = 1e-16\n",
    "\n",
    "A[:,1] = A[:,0] + delta*A[:,1]\n",
    "print(\"cond = %g\" % np.linalg.cond(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a problem we know the answer to:\n",
    "\n",
    "Let $x={\\bf 1}$, then $x$ solves the problem\n",
    "$$\n",
    "A x = b\n",
    "$$\n",
    "where $b = A {\\bf 1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This is the exact solution\n",
    "xexact = np.ones((n,))\n",
    "b = A.dot(xexact)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This is the approximated solution\n",
    "xnum = np.linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we are solving with LU with partial pivoting, the residual should be small!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Residual Versus Error\n",
    "$$\n",
    "r = b - A x\n",
    "$$\n",
    "whereas\n",
    "$$\n",
    "e = x_{exact} - x\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xexact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xnum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r = b - A@xnum\n",
    "e = xexact - xnum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"norm of residual = \", la.norm(r))\n",
    "print(\"norm of the error = \", la.norm(e))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The condition number of A is high (ill-conditioned problem), and hence the error bound is also high."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
